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Abstract 

In this work we propose a generalization of the Moment Guided Monte Carlo method 
developed in [11]. This approach permits to reduce the variance of the particle methods 
through a matching with a set of suitable macroscopic moment equations. In order to guarantee 
that the moment equations provide the correct solutions, they are coupled to the kinetic 
equation through a non equilibrium term. Here, at the contrary to the previous work in 
which we considered the simplified BGK operator, we deal with the full Boltzmann operator. 
Moreover, we introduce an hybrid setting which permits to entirely remove the resolution of the 
kinetic equation in the limit of infinite number of collisions and to consider only the solution 
of the compressible Euler equation. This modification additionally reduce the statistical error 
with respect to our previous work and permits to perform simulations of non equilibrium gases 
using only a few number of particles. We show at the end of the paper several numerical tests 
which prove the efficiency and the low level of numerical noise of the method. 

Keywords: Monte Carlo methods, hybrid methods, variance reduction, Boltzmann equation, 
moments closure, fluid equations. 



1 Introduction 

The kinetic description of a gas is based on the Boltzmann equation [21 [5] , which is a linear hyper- 
bolic partial differential equation with a source term. This equation is satisfied by the probability 
distribution function of the particles depending on the time t > 0, the space variable x € M. d and 
the velocity v € R d . In addition to the high dimensionality of the problem, the Boltzmann colli- 
sion term that characterizes the kinetic equation is a five fold integral very hard to treat due to its 
nonlinear nature and to the physical properties which need to be preserved during its resolution. 

For these reasons, the numerical simulation of the Boltzmann equation with deterministic 
techniques appears to be very difficult and thus, in practice, probabilistic techniques such as 
Direct Simulation Monte Carlo (DSMC) methods are extensively used in real situations due to 
their large flexibility and low computational cost compared to finite volume, finite difference or 
spectral methods for kinetic equations [TJ [21 El El HEl US] • On the other hand, DSMC solutions 
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are affected by large fluctuations. Moreover, in non stationary situations the impossibility to use 
time averages to reduce fluctuations leads to, either poorly accurate solutions, or computationally 
expensive simulations. To overcome this problem, several methods have been developed. We quote 
[B] for an overview on efficient and low variance Monte Carlo methods. For applications of variance 
reduction techniques to kinetic equation we mention the works of Homolle and Hadjiconstantinou 
[22] and [23] . We mention also the research papers by Pullin [30] which developed a low diffusion 
particle method for simulating compressible inviscid flows and the work of Weinan E. and co- 
authors [20] on the development of multiscale numerical methods. We quote also the work of 
Lcmou and co-authors who developed a low variance method for the Vlasov equation close to the 
fluid limit [9] . We finally recall the results obtained by the author of the present paper [15l [16] on 
the construction of efficient and low variance methods for kinetic equations in transitional regimes. 

A second important drawback of kinetic approaches is that the collision term becomes stiff 
close to the fluid regimes. A non dimensional measure of the importance of collision is given by 
the Knudsen number which is large in the rarefied regions and small in the fluid ones. Thus, 
standard computational approaches lose their efficiency due to the necessity of using very small 
time steps in deterministic schemes or, equivalently, a large number of collisions in probabilistic 
approaches. One possibility which permits to avoid the severe time step restrictions caused by the 
small free path limit is to construct the so-called Asymptotic Preserving schemes [TBI El] which 
permits to avoid time steps dependence of the mean free path. Another possibility to avoid the 
excessive computational cost and the stiff regimes is to construct numerical schemes which combine 
continuum models with microscopic kinetic models [H [5] [lOj [H] [13j [M] [21] [26] [31] [32]. However, 
this approach has several drawbacks, first the identification of the different regions is not trivial 
and second match the two models at the interfaces is not simple. For these reasons alternatively 
approaches are often preferred. 

The two main goals of this work are first, design a particle closure of the macroscopic fluid 
equations which avoid the severe time step restrictions caused by the collisional scale and second 
to reduce the statistical error caused by the use of DSMC methods. This approach has explained 
in details in [13] can also be the basis of very efficient domain decomposition strategies. In our 
first paper [11], we developed a method for the simple case of the BGK equation which consists 
in forcing particles to match prescribed sets of moments given by the solution of deterministic 
equations. Now, we generalize the approach to the case of the Boltzmann equation. Moreover, 
we construct the scheme in such a way that the time step restriction due to the small Knudsen 
number are avoided, thanks to the introduction of exponential schemes for the resolution of the 
Boltzmann collision term [17]. In particular, in the limit e — > the scheme becomes automatically 
an high order method for the compressible Euler equation in which the numerical noise completely 
disappears. In this sense, the method belongs to the class of the so-called asymptotic preserving 
methods [13 [18j UJ [24] . 

The method is based on the following idea. The resolution of the problem is done through the 
kinetic equation which is solved by DSMC method and by a set of closed moment equations which 
are solved by means of finite volume techniques. In order to provide the correct solution for all the 
regimes of the Knudsen number, the moment equations are coupled to the Boltzmann equation 
through a kinetic correction term, which takes into account departures from thermodynamical 
equilibrium. In this sense, we perform a closure to the infinite set of moment equations, using 
the solution of the kinetic equation. Observe that the kinetic model will be sufficient to give the 
correct solution to the problem. However, as already explained, the Boltzmann equation should 
be resolved by very fast methods because of its too high computational cost. Typically in real 
situations Monte Carlo is the only possible choice. The drawback is that the solution is only 
poorly accurate. On the other hand, the moment equations can be resolved by using high order 
finite volume techniques. The result is a method which has a computational cost comparable to the 
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cost of DSMC method but which a smaller level of numerical noise. Moreover, the stiffness of the 
problem is removed thanks to the exponential resolution of the collision term. This gives a method 
which is faster than classical DSMC methods [TJ [5] in the limit of the Knudsen number which goes 
to zero i.e. in the fluid limit. Finally, in order to recover at each time step the same moments 
from the solution of the kinetic equation and from the solution of the moments equations and then 
advance in time, we constrain the DSMC method to match the moments obtained through the 
deterministic resolution of the macroscopic equations in such a way that the higher accuracy of the 
moments resolution improves the accuracy of the DSMC method. We experimentally show that 
this is indeed the case. 

The remainder of the paper is organized as follows. In the next section we recall some basic 
notions on the Boltzmann equations and its fluid limit. The scheme is described in section 3. 
In particular, the numerical method used for the kinetic equation is described in section 3.1, the 
numerical method used for the moment equation is described in section 3.2, while the matching 
procedure is described in section 3.3. In section 4 numerical examples which demonstrate the 
capability of the method are presented. Finally some conclusions and remarks are drawn in the 
last section. 



2 The Boltzmann equation and its fluid limit 

We consider the Boltzmann equation of rarefied gas dynamics [8, 

d t f + v-V x f = Q(fJ) (1) 

with initial data 

f\t=o = fo- (2) 

Here f(x, v, t) is a non negative function describing the time evolution of the distribution of particles 
with velocity »el 3 and position x <G fl C at time t > 0. In the sequel for notation simplicity 
we will omit the dependence of / from the independent variables x,v,t unless strictly necessary. 
The operator Q(f, f) which describes particle interactions has the form 

Q(fJ)= [ B(\v-v*\,n)[f(v')f(vl)-f(v)f(v*)]dv*dn (3) 

JR 3 xS 2 

where 

v 1 =v + -(v- u*) + -\v - u*|n, = v + -(v - v*) - -\v — v*\n, (4) 

and B(\v — u*|,n) is a nonnegative collision kernel characterizing the microscopic details of the 
collision given by 

B(\v-v*\ ) n) = *(^^ r n)\v-v*\i, (5) 

with 7 6 [0, 3) and a the collision cross section which in general depends on the collision angle. The 
case 7 = 1 is referred to as hard spheres case, whereas the simplified situation 7 = 0, is referred to 
as Maxwell case. Note that in most applications the angle dependence is ignored and a is assumed 
constant [2]. 

The operator Q(f, /) is such that the the local conservation properties are satisfied 

mQ(f,f)dv=:(mQ(fJ)) = (6) 
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where m(v) = (l 1 v 1 -Sp) are the collision invariants. In addition it satisfies the entropy inequality 

^/ f log fdv= f Q(fj)\ogfdv<0. (7) 

Integrating |lj against its invariants in the velocity space leads to the following set of non closed 
conservations laws 

d t (mf)+V x -(vmf)=0. (8) 

Equilibrium functions for the operator Q(f, /) (i.e. solutions of Q(f, /) = 0) are local Maxwellian 
of the form 

P f — \u~v\ 2s 



where p, u, T are the density, mean velocity and temperature of the gas in the x-position and at 
time t defined as 

(p,pu,pe) = (mf), E = pe, T = -^-(E ~ p\u\ 2 ). (10) 

We will denote by 

U=(p,pu,E), M[U] = M f , (11) 

clearly we have 

U={mM[U]}. (12) 

Now, when the mean free path between particles is very small the collision rate is large and we 
can rescale the space and time variables in ((T|) as 

x = ex, t = et (13) 

to obtain 

d t f + v-V x f=±Q(fJ) (14) 

where e is a small parameter proportional to the mean free path, the so-called Knudsen number, 
and the primes have been omitted to keep notation simple. 

Passing to the limit for e — > leads to / — > M[U] and thus we obtain the closed hyperbolic 
system of compressible Euler equations for the macroscopic variables U 

8 t U + V x -F(U)=0 (15) 

with 

F(U) = (vmM[U]) = {pu, gu(g> u + pi, Eu + pu), p = pT, 
where / is the identity matrix. 

3 The Monte Carlo method for the Boltzmann equation 

The starting point of the method is the following micro-macro decomposition 

f = M[U)+g. (16) 

The function g represents the non-equilibrium part of the distribution function. From the definition 
above, it follows that g is in general non positive. Moreover since / and M[t/] share the same first 
three moments we have 

(mg) = 0. (17) 
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Now U and g satisfy the coupled system of equations 



d t U + V x -F(U) + V x -(vmg) 



(18) 
(19) 



dtf + v- V x f 



e 



We skip the proof of the above statement and refer to [12] for details on the decomposition of the 
distribution function and the different coupled systems which it is possible to derive. 

As explained in the introduction, the goal of the method is to solve equation (fT5|) , which means 
to close the system of moments equations. In order to accomplish that task, we need to know the 
time evolution of g, which is given by the solution of equation (|19[) . Of course solving equation 
(TT9"]) will be sufficient to compute the solution also of (ITS)) . However, the numerical solution of 
equation (p~9|) is very expensive and thus we seek for very fast solvers, as for instance the Monte 
Carlo method. The counterpart of fast solvers is that the solution of (fT9"| will be normally known 
with very low accuracy. Thus, our statement is the following: if the departure from equilibrium 
g is small, the low accuracy and the high level of numerical noise with which g is known from 
(|19[) will still permit to compute the moments in (I18|) with an error which will be smaller than the 
error with which equation (fl9|) is solved. We will numerically show that this is indeed the case. In 
particular, thanks to the way in which the collisional operator is solved, the coupling method will 
permit to make disappear g from the computation of the moments in the limit e —> 0. The precise 
time marching procedure will be detailed next, the method can be summarized as following: at 
each time step t n 

1. Solve the kinetic equation (fT9"|) with a Monte Carlo scheme and obtain a first set of moments 



2. Solve the fluid equation (fT5|) with a finite volume/difference scheme using particles to close the 
system and obtain a second set of moments U n+ . This means evaluate the term X7 X ■ (vmg) 
with particles. 

3. Match the moments of the kinetic solution with the fluid solution through a transformation 
of the samples values f n+1 = T(f*) so that (mf n+1 ) = U n+1 and continue. 

In principle, for Step 1, one can use any Monte Carlo method or more generally any low accurate 
but very fast solver. Step 2 and 3 of the above procedure require great care since they are the 
key points of the method, they involve the evaluation of Va; • (vmg) and the moment matching 
procedure. 

Finally let us note that, in principle, it is possible to improve the method, adding to system 
(|18[) additional equations for the time evolution of higher order moments and get 



with 77i„ = v n and n > 3. However, we will not discuss this possibility here. In the sequel we give 
the details of the above three Steps. 

3.1 Solution of Boltzmann equation by exponential Runge-Kutta Monte 
Carlo methods 

The starting point in the solution of kinetic equations by Monte Carlo methods is given by an 
operator splitting of (TTJ) or equivalently (|T9l in a time interval [0, At] between relaxation 



U* = (m/*>. 



dt(m n f) + V x ■ (vm n f) = {m n Q{f, /)) 



(20) 



d t f = -Q(f,f), 



(21) 
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and free transport 



dtf + v- V x f = 0. 



(22) 



Even if this splitting is limited to first order it permits to treat separately the hyperbolic free 
transport from the stiff relaxation step. Observe that high order splitting schemes are possible 
with Monte Carlo techniques. However, we did not consider this possibility in this work, we 
remind to the future for the exploration of high order time discretizations. 

We introduce now a space discretization of mesh size Air and a time discretization of time step 
At. The discretization of the domain is not needed for the transport step which is solved exactly 
by pushing the particles. However, it is necessary to solve the collision part of the problem which 
acts locally in space and thus it turns to be necessary to solve the full problem. In Monte Carlo 
simulations the distribution function / is discretized by a finite set of N particles 

N 

/ = ^ E 6 ( x - Mt))S(v - Vi{t)), (23) 

i=l 

where Xi(t) represents the particle position in the three spatial directions, Vi(t) the particle veloc- 
ities in the velocity space, and cti(t) the weight to associate to each particle. The constant m p is 
defined at the beginning of the computation in the following way 

m p = i I I /(* = °1 dvdx - ( 24 ) 

During the transport step (|22|) . the particles move to their next positions according to 

Xi(t + At) = Xi(t)+Vi(t)At (25) 

where (|23|) with ([25]) and a constant oti{t) represents an exact solution of equation (|22]l. The role 
of the function a will be clarified in the moment matching procedure. 

We consider now the solution of the collision step (|2T|) . The effect of this step is to change the 
shape of the velocity distribution and to project it towards the equilibrium distribution leaving 
unchanged mass, momentum and energy. The collision operator acts locally in space which means 
that we solve it independently in each spatial cell. The particles are assumed to have all the same 
weight in one cell. The initial data of this step is given by the solution of the transport step because 
of the splitting. 

There exists many different way to solve the collisions, see for instance [TJ|5J|S], however, we seek 
for a method which will be stable for all choices of At independently of e. This will permit to avoid 
the stiffness of the equation and at the same time to reduce the statistical fluctuations. Since we 
aim at developing unconditionally stable schemes, the most natural choice would be to use implicit 
solvers applied to (|2"Tj) . Unfortunately the use of fully implicit schemes for (|2"Tj) is unpracticable 
in the case of DSMC methods. To overcome these difficulties, we introduce a method based on 
exponential integrators. First we rewrite the homogeneous equation (|21[) in the form 

d t f=±(P(f,f)- f if), (26) 

where P(f, f) — Q(f, f) + and [i > is a constant such that P(f, f) > 0. Typically fi is 
an estimate of the largest rate of the negative term in the Boltzmann operator. For example for 
Maxwellian molecules we have 

P(f,f) = Q + (f,f)(v)= [ [ b (cos6)f(v')f(vl)dojdv*, (27) 

JS 2 



G 



and n = q. By construction the following property holds 



-(mP(fJ)) = (mf) = U, (28) 
M 

this means that P(f, f)/n is a density function. The homogeneous equation can be written now, 
adding and subtracting M[U], in the form 

d t f = £ {^P- - M[U]\ + ^{M[U\ - /). (29) 

Note that even if M[U] is nonlinear in / thanks to the conservation properties ©, it remains 
constant during the relaxation process. 

We apply now to the reformulated equation ([29]) an approach based on the so-called exponential 
integrators where the exact solution of the linear part is used for the construction of the numerical 
method [17] ■ In order to derive the methods it is useful to rewrite (|2"91 as 

at e 

The above form is readily obtained if one multiplies ([29| by the integrating factor exp(fit/e) 
and takes into account the fact that M[J7] does not depend of time during the collisional process. 
A class of explicit exponential Runge-Kutta schemes is then obtained by direct application of 
an explicit Runge-Kutta method to (|30p . More generally we can consider the family of methods 
characterized by 



3=1 

l- e -CiHAt/e\ M pn^ • = ^ 



(31) 



pn+1 



1 - e ~^ At/£ ) M[U n 



(32) 



where /* is the distribution function value after the transport step, are called stages, Ci > 0, 
while the coefficients Ay and the weights Wi are such that 

Aij(0) = a,ij, Wi(0) = Wi, i,j = l,...,v 

with coefficients aij , a and weights Wi given by a standard explicit Runge-Kutta method called the 
underlying method. Various schemes come from the different choices of the underlying method. 
The most popular approach is the integrating factor (IF) method. For the so-called Integrating 
Factor methods, which correspond to a direct application of the underlying method to (|30l) . we 
have 

Aij{\) = a, J e" (c ^ c,)A , i,j = l,...,v, i>j 

(33) 

Wi(X) = Wie-^- c ^ x , i = l,...,v, 
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with A = [lAt/s. In the sequel, among all possible choices for discretizing (|2"6")1 . we will use the 
first order in time integrating factor scheme which reads 



r+1 = r + ^ e _^ ^pif^n _ M[lfn] \ + Mn (34) 

which is based on the simple firs order in time explicit Euler Runge-Kutta method. Observe that 
this method permits to avoid the stiffness of the collisional operator. This is, in fact, unconditionally 
stable for all choices of time step At. At the same time, the method proposed is, in the limit e — > 0, 
nothing else than a kinetic scheme for the numerical solution of the compressible Euler equation. 
Observe, in fact, that when e — > at each time step the distribution function is projected on the 
equilibrium distribution M[U\. 

The Monte Carlo interpretation of equation (f34|) is the following. If we define 

uAt a At 11 At „ / uAt gAt 11 At \ „ 

A = e~ — , B=!—e~—, C = ( 1 - [ — e -— _ e ~— j ) C = 1 — A — B, (35) 

with probability A the velocity of the particle does not change during the collision step. With 
probability B the particle occurs in one collision, the details of the collision are described by the 
operator P, this is, for instance, the gain part of the collisional process, described by Bird [2]. With 
probability C the velocity of the particle is replaced by a particle sampled from the Maxwellian 
distribution M [U]. We will discuss in the matching section the way in which we link the solution 
of the kinetic equation to the solution of the moment equations. 



3.2 Solution of the moments equations 

In this section we discuss the discretization of the moments equations. We consider, in order 
to write the discretization, the function g n known at time n. Our scope, in the construction of 
the numerical scheme is to take advantage from the knowledge of the Euler part of the moment 
equations 

d t U + V x -F(U)+V x -{vmg) = 0. (36) 

S v ' 

Euler equations 

Thus, the method is based on solving first the set of compressible Euler equations, for which many 
efficient numerical methods already exsist, and then considering the discretization of the kinetic 
flux V x • (vmg). To that aim, for the space discretization of the compressible Euler equations we 
use a second order MUSCL central scheme, while a backward discretization is used for the time 
derivative in all cases. For simplicity, we indicate in the same way the numerical flux in one or in 
more spatial dimensions. Thus, we have 



u;-uf t Vi+i/2(tf")-^-i/2(tf") 



At 

The discrete flux reads [25l [27] 

l / i / -m- -r-ft \ • • / -r -r-tt \ \ ^ / ~r -f <Y1 -r tVI \ 



0. (37) 



^•+i/2(t/") - ^(F(U?) + F(U? +1 )) -a{U? +1 U?) + -(a^ (38) 



where 



= {F{U- +1 ) ± aUf +l F{UV<) T aUf) y{x^) (39) 
with if the Van Leer slope limiter 

<P(X) = T~' (40) 
1 + X 
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finally, the variable \ ± is defined as following 



± = F(U?)± a U?-F{U?_ 1 ) T *U?- 1 

Xj T?tTTn. \ J TTn T?(TTn\ t~ „TTn \ I 



F{Uf +1 ) ± aU? +1 - F(U}) T all? 



where the above ratio of vectors is defined componentwise and a is equal to the largest eigenvalue 
of the Euler system. 

We now discuss how to discretize the non equilibrium term V x - < vmg >. To this aim, the 
first step is to introduce a filter to eliminate some fluctuations. We choose the weighted moving 
average method, which is a convolution of the pointwise value of g with a fixed weighting function. 
Thus, given the function g n , we define 

j k=K k=K 

& = 2K + 1 £ w ^- fc ' E Uk = L ( 42 ) 

k=-K k=-K 

The smoothing permits to reduce the fluctuations but on the other hand it causes a degradation of 
the accuracy in the solution. Thus, in practice, we use a weak filter in order to keep the solution as 
precise as possible. In the numerical simulations, we used K = 1, w_i = lo\ = 1/6 and wo = 2/3. 
Observe that a different smoothing can be used at the particle level instead that at the moments 
level. This different smoothing has to be applied during the reconstruction of the moments of g 
from the particles. We discuss this possibility in the next section. Once that the filer is applied, the 
non equilibrium term is discretized with the same second order MUSCL scheme used for computing 
the flux of the Euler equations where, however, the numerical diffusion is taken equal to zero. This 
is because this term is at the leading order a diffusion term, and thus it does not need numerical 
diffusion for stability to be assured. Thus the final scheme, for the moment equations, reads 

Vj-Vj | M1^E) =0 , (43) 



At Ax 

U ? +1 U _j ^+i/2(< vmr >) - ^--i/ 2 (< vmr » = 
At Ax ' 

where ^j+i/2(< vmg n >) is defined as 

*j+i/ 2 (< vmg n >) = i(< vmg 1 } > + < wm<# +1 >) + ^(a^ - <7™ ij+1 ) (45) 

with 

a g,j = (< K+i > - < vm 9? >) ( 46 ) 

and 

x n = < vmgj > < vmg^ > 
9 ' J < vmg™ +1 > — < vmg} > 

where again the above ratio of vectors is defined componentwise. 
3.3 The Moment Matching 

In this section we discuss the details of the moment matching and the coupling between the kinetic 
and the macroscopic model. Suppose, known / n_1 , <?™ _1 and U 11 ^ 1 . The time marching procedure 
is the following 
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1. 



Solve the moments equations from time n — 1 to time n using equations (|43fl44|) this gives 



Solve the transport part of the kinetic equation from time n — 1 to the intermediate stage * 
pushing the particles as in equation (|25[) . 



3. Compute the moments of the Boltzmann equation after the transport of the particles. This 
gives U n , the collision part, being conservative, does not alter the moments. The reconstruc- 
tion of the moments from the particles will be discussed next. 

4. Match the moments of the kinetic and the macroscopic equations forcing the particles to 
have the moments U n , i.e. U n —> U n . 

5. Compute the collision step of the kinetic equation through ([34]) this gives f n . 

6. Use the computed value of /" to compute g n . 

7. Compute g n with the moving average technique (equation I42[) . 

8. Plug the value of g n in the moments equations and compute U n+1 using (|43M4j) and continue. 

We discuss now, how to reconstruct the moments from the particles (point 3 of the above 
procedure), the matching of the different moments (point 4) and the computations of g n from /" 
(point 6). 

We first consider the matching of the mass. To this aim, let consider the set of particles 
X\, . . . , Xn x . inside the cell Xj after the transport step. The corresponding mass is computed as 

N x - 

0j=J2m p ar\ J = l,N x (48) 

others possible reconstructions will be discussed next. In our previous paper, among the 

possible techniques that can be used to restore a prescribed density we choose to replicate or 
discard particles inside the cells. Here, in order to restore the mass we assign to each particle 
inside the cell j the new value 

p" 

a? = jf-Vi = l,..,N Xj , (49) 

where g" is density computed from the solution of the moments equations. This rescaling of the 
mass introduce an error in the evaluation of the distribution function /. However, we recall that the 
two models, the microscopic and macroscopic one, give the same solution in term of the moments 
apart from the numerical errors. In particular, the DSMC method gives solutions which oscillates 
around the exact solution of the problem. This implies that, the above renormalization of the mass 
only force the density to be closer to the value furnished by the exact solution of the problem. In 
fact, the moments equations furnish solutions which contain less fluctuations with respect to the 
DSMC method. The reason for this lower level of noise of the macroscopic equations is that only 
the perturbation from the equilibrium g is computed statistically and not the full solution as in 
the original DSMC method. 

We discuss now the matching of the momentum and of the energy, which are done after the 
matching of the density. In the Monte Carlo setting these moments can be obtained by the following 
piecewise constant reconstructions 

N Xj N X] 

ui=—k^y^Vi ^ ^—^y^?, (so) 
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where Uj is a vector representing the mean velocity in the three spatial directions. To match mean 
velocity and energy with those of equations (|43H44I) we then apply the transformation described in 
[B] which permits to get a new set of velocities V* defined by 



(u n ) 2 

V* = (V i -u«)/c + u? c=J %_^ 2 , i=l,...,N Tj (51) 



which gives 



N Xj N Xj 

Vk* = u?, -J— V hv;\ 2 =e r - 



After the matching at time n, the next point in the time marching procedure described before, 
regards the collisions. This step has been already described in section [3TTI It uses equation (|34| 
to compute from /*, which is now the distribution function value after the transport and the 
matching, the new distribution f n . 

We now discuss the sixth point of the method: the computation of g n from /". The perturbation 
g n is given by 



jn-i + ^ er ^^l Li L + 1 - e-~ - t—e-^- M\m] - M P" 

e // \ e 



= r- + f^ e -^nr-\r- 1 ) _ (-^ + /^ e -^\ Mm (52) 

£ fl V £ J 

where once again we used the first order in time integrating factor scheme (|34[) . The above 
expression tells us that the moments of <?" can be obtained as a contribution of two terms. The 
first term is obtained by reconstructing from the particles the moments as in (|48l [50]) . The second 
term is obtained by integrating over the velocity space the analytic expression of the Maxwellian 
distribution. This implies directly that the moments related to the second part does not contain 
any statistical error. Finally, observe that in the limit e — > the contribution of the perturbation 
g goes to zero. As a consequence, the method proposed resolves a macroscopic system which is 
as closer as the Knudsen number diminishes to the compressible Euler equation. In the limit of 
£ = 0, the scheme exactly solve the compressible Euler equations without any source of statistical 
error. Finally, the discrctizcd moments equations (I44[) can be rewritten as 

U; +1 - U; = gj+i/a(< vm 9 n >) - *j-i/2(< vrng" >) = 
At Ax 

*W« vm (e-^ + i^e"^) Q - M)« » 



Ax 

%-i/ 2 (< vm (e"^ + BM e -^ Q _ M T >) 



Aa 



(53) 



where 



e -zM , n _i + ^At e -afi pjrz\r^) 

<« = (psr^Ff (54) 

As Ai/e grows, which means that the system approaches the equilibrium, the contribution of the 
kinetic term vanishes even though it is evaluated through particles. Observe that, this does not 
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happen if we just compute the kinetic term • (vmg) from the particles without considering 
the structure of the distribution function /. This dramatically decreases fluctuations when the 
Knudsen number is small. 

As a conclusion for this section, we discuss some different reconstructions of the moments 
starting from the particles. Instead of using the piecewise constant reconstructions (|48fl50[) . one 
can think to smoother reconstructions which can probably further diminish the statistical noise. 
This procedure will be alternative to the moving average algorithm used for computing g from g 
(142p . In other words, we first smooth the contribution of the particles in the computation of the 
macroscopic quantities and we then after just apply the finite volume method to set of resulting 
moments equations without using average algorithms. 

A general way to compute the integral over the velocity space is to sum over the particles with 
the use of weight functions. In this strategy, the value of the perturbation g will be given by 

N 

(vmg) = ^2 m p a i v i B(X l - Xj)rm j = 1, .., N x (55) 

i 

where B > is a suitable weight function s.t. 

/ B(x)dx = 1. 
Jr 

For example, B(x) = 1 if |x| < Ax/2 and B(x) = elsewhere, gives rise to the previous piecewise 
constant reconstruction, while B{x) = 1/ Ax if |a;| < Ax/2 and B(x) — elsewhere, is the so called 
Nearest Grid Point procedure in plasma physics [3] . Smoother reconstructions can be recovered by 
convolving the samples with a bell-shaped weight like a B-spline. Note that the value B(x) has a 
strong influence on the fluctuations in the reconstructed function, and in general should be selected 
as a good compromise between fluctuations and resolution as for the moving average method. 

4 Numerical results 

In this section we report some numerical results for the Moment Guided method on different test 
cases. First, we consider two classical shock problems and then, we perform an accuracy test using 
a smooth periodic solution . For the two shock tests, we compare the moment guided (MG) solution 
with a Monte Carlo (MC) solution which employs the first order exponential Runge-Kutta method 
together with a splitting technique between the transport and the collision parts. We report also 
in each figure the results for the limit compressible Euler equation using the same second order 
MUSCL scheme which has been used for the equilibrium part of the Moment Guided method. 
Finally, the reference solutions for the two shock test problems are obtained by using the same 
Monte Carlo method described above with the same number of cells, in which, however, the number 
of particles is such that the statistical noise is very small. 

All the tests considered are relative to unsteady problems. This choice reflects the fact that the 
method proposed in this paper is specifically aimed at situations in which the classical variance 
reduction techniques which employs time averaging, typical of DSMC methods, cannot be used or 
turns out to be useless, since time- aver aging or using more particles leads to the same computational 
effort. 

Unsteady shock test 

The first problem concerns the study of an unsteady shock. The figures [1] to [3] consider the same 
initial data for the density g = 1, the mean velocity u = — 1 and the temperature T = 1 for 



12 



Density, Knudsen number=0.01 Density, Knudsen numbered. 01 




0.5 1 1.5 0.5 1 1.5 



Velocity, Knudsen number=0.01 Velocity, Knudsen number=0.01 




Temperature, Knudsen number=0.01 Temperature, Knudsen number=0.01 




Figure 1: Unsteady shock test: Solution at t = 0.18 for the density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MG method (right). Knudsen number 
e = 10 -2 . Reference solution: dash dotted line. Euler solution: continuous line. Monte Carlo or 
Moment Guided: circles plus continuous line. 400 particles per cell. 
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Density, Knudsen number=0.001 Density, Knudsen number=0.001 




0.5 1 1.5 0.5 1 1.5 

X X 



Velocity, Knudsen number=0.001 Velocity, Knudsen number=0.001 




Temperature, Knudsen number=0.001 Temperature, Knudsen number=0.001 




Figure 2: Unsteady shock test: Solution at t = 0.18 for the density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MG method (right). Knudsen number 
e = 10 -3 . Reference solution: dash dotted line. Euler solution: continuous line. Monte Carlo or 
Moment Guided: circles plus continuous line. 400 particles per cell. 
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Density, Knudsen number=0.0001 Density, Knudsen number=0.0001 




0.5 1 1.5 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

X X 



Velocity, Knudsen number=0.00u1 Velocity, Knudsen number=0.0001 




Temperature, Knudsen number=0.0001 Temperature, Knudsen number=0.0001 




Figure 3: Unsteady shock test: Solution at t = 0.18 for the density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MG method (right). Knudsen number 
e = 10 -4 . Reference solution: dash dotted line. Euler solution: continuous line. Monte Carlo or 
Moment Guided: circles plus continuous line. 400 particles per cell. 
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different initial Knudsen number values, ranging from e = 10~ 2 to e = 10~ 4 . Specular reflection 
boundary conditions are used on the left side to produce the shock wave. The right boundary 
condition prescribes inlet flow. The number of cells is 150 while the time step is given by the 
minimum of the ratio of Ax over the maximum velocity owned by the particles and the ratio of Ax 
over the larger eigenvalues of the compressible Euler system. The final time is 0.18. The Knudsen 
number value does not play any role in the choice of the time step, being the method asymptotic 
preserving and thus independent from the small scale constraint prescribed by e. Each figure 
depicts the density, the mean velocity and the temperature from top to bottom, with the pure 
Monte Carlo solver on the left and the Moment Guided method on the right. At the beginning of 
the simulation 400 particles per cell are used for both methods Monte Carlo and Moment Guided, 
the solutions intentionally still contain some fluctuations. This is done to clearly show the difference 
in term of statistical error between a Monte Carlo method and our method. In addition, we report 
the solutions of the compressible Euler equations and the reference solution computed by the same 
Monte Carlo method where 5 10 5 particles per cell are employed. 

The three figures [l|[3] show a large reduction of fluctuations for all cases analyzed. Observe 
that the reduction of the fluctuations depends on the Knudsen number value. In particular in 
the case e = 10 -4 the solution is very close to the limit solution. Observe also that the solution 
furnished by the Moment Guided method when e = 10~ 4 is closer to the shock with respect to 
the reference Monte Carlo solution. This can be interpreted as an over relaxation problem of the 
Moment Guided method. However, an in deep simulation analysis we did not report, shows that, 
the difference between the MC and the MG methods when e = 10~ 4 is due to the more diffusive 
behavior of the Monte Carlo method. This is caused by the large time steps At allowed by the 
exponential method. These large time steps introduce a numerical diffusion in the schemes which 
is more important for the Monte Carlo method and less important for the Moment Guided one. 

Sod shock tube 

The second problem analyzed is the Sod tube test. Again this choice relies on the fact that the 
method is specifically aimed to the study of unsteady problems. In this case, we want to study 
the capability of the method in describing different waves with respect to the simple shock wave. 
The figures 0] to [7] consider the same initial data for the density g = 1 upstream and g — 0.125 
downstream of the initial shock, the mean velocity u = in all the domain and the temperature: 
T = 5 upstream and T = 4 downstream of the shock. Different initial Knudsen number values are 
considered in the test, which range from e = 10~ 2 to e = 10~ 4 while the number of cell is 200. 
In figures [U [6] and [7] 200 particles per cell are also used on average, more precisely 200 particles 
are present at the beginning of the simulation and then distributed accordingly to the density in 
each cell. The final time is Tfi n = 0.08, while the time step as previously is given by the minimum 
between the the ratio of Ax over the maximum velocity of the particles and the ratio of Ax over 
the larger eigenvalue of the compressible Euler equations. Again the Knudsen number value does 
not constrain the time step, thanks to the exponential method employed in the resolution of the 
collision integral. 

As for the unsteady shock test each figure depicts the density, the mean velocity and the 
temperature from top to bottom, with the pure Monte Carlo solver (left) and the Moment Guided 
method (right). The reference solution is obtained through the Monte Carlo method in which the 
number of particles is taken equal to 2 10 7 with a solution averaged 5 times, this gives a final 
number of 10 s particles employed for constructing the reference solution. The figures show good 
results for all ranges of Knudsen numbers in terms of reduction of fluctuations. In particular, we 
clearly see a strong reduction of fluctuations for e — 10~ 3 (figure [g}, while for e = 10~ 4 (figure E} 
the fluctuations are almost completely disappeared. In the case of e = 10 -2 (figure 3}, we see that 
even if the statistical error is smaller for the Moment Guided method than for the Monte Carlo 
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Temperature, Knudsen number=0.01 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.( 



Figure 4: Sod shock tube test: Solution at t = 0.08 for the density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MG method (right). Knudsen number 
e = 10 -2 . Reference solution: dash dotted line. Euler solution: continuous line. Monte Carlo or 
Moment Guided: circles plus continuous line. 200 particles for cell. 
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X X 



Figure 5: Sod shock tube test: Solution at t — 0.05 for the density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MG method (right). Knudsen number 
e = 10 -2 . Reference solution: dash dotted line. Euler solution: continuous line. Monte Carlo or 
Moment Guided: circles plus continuous line. 1000 particles for cell. 
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Density, Knudsen number=0.001 Density, Knudsen number=0.001 
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Figure 6: Sod shock tube test: Solution at t = 0.08 for the density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MG method (right). Knudsen number 
e = 10 -3 . Reference solution: dash dotted line. Euler solution: continuous line. Monte Carlo or 
Moment Guided: circles plus continuous line. 200 particles for cell. 
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Figure 7: Sod shock tube test: Solution at t = 0.08 for the density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MG method (right). Knudsen number 
e = 10 -4 . Reference solution: dash dotted line. Euler solution: continuous line. Monte Carlo or 
Moment Guided: circles plus continuous line. 200 particles for cell. 
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method, both solutions are still far from the converged solution indicated by the dash-dotted line. 
Thus, we reported the same simulation in figure [5] in which 1000 particle per cell on average are 
employed for both methods. This figure permits to show the faster convergence of the Moment 
Guided method towards the reference solution also in the case of larger Knudsen. Finally, thanks 
to the high order discretization of the moments equations, observe that when the gas is close to 
the fluid limit (figure [7]), the shock is very well represented. 



Accuracy test 

In this final part, we report the results of a stochastic error analysis with respect to the number 
of particles. Here, as reference solution we considered the average of M independent realizations 



U 



MC 



1 M 

TjE^M/c (56) 

k=l 



and 



- 1 M 

Umg = JjJ2 ^M/g (57) 
fe=i 

where the two subscripts MC and MG indicate respectively the reference solution for the Monte 
Carlo method and for the Moment Guided method. We used two different reference solutions since 
the two schemes present different discretization errors and thus they converge, when the number 
of particles goes to infinity, to different discretized solutions. The difference between the two limit 
solutions (Monte Carlo and Moment Guided Monte Carlo) is mainly due to the different numerical 
diffusion introduced by the two methods. This has been made clear for instance in the unsteady 
shock test case when e = 10 -4 . In this already explained, the Monte Carlo solution is 

more diffusive with respect to the Moment Guided one and thus, in order to measure only the 
statistical error, we cannot consider the same reference solution for the two schemes. The two 
reference solutions are obtained by fixing the time step and mesh size and letting the number of 
particles goes to infinity. The number of samples used to compute the reference solution is on 
average 10 5 per cell, while the number of realizations is M=5. The number of spatial cells is 100, 
the time step is fixed to 10 -3 for all the values of the Knudsen number. Taking a large number of 
samples for computing the reference solutions gives negligible stochastic errors. However, at the 
same time, both solutions involve space and time discretization errors. In any case, the amount of 
such errors does not change when the number of particles varies. Therefore by comparing solutions 
obtained with a given At, Ax, but with a finite number of particles to reference solutions obtained 
with the same At, Ax, but with a very large number of particles, we obtain a true measure of the 
error originating from the stochastic nature of the method. 
In practice we measured the quantity 



M N 

s2 w = mEE(^-^) 2 ( 58 ) 

where Uj represents the reference solutions, N — 100 and M — 10. The test consists of the 
following initial data 

2jtX 

g(x, 0) = 1 + a e sin — — 

^TTX 

u(x, 0) = 1.5 + a u sin — — (59) 

Li 
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Figure 8: Statistical error test: Solution at t = 0.05 for density (top), velocity (middle) and 
temperature (bottom). MC method (left), Moment Guided MC method (right). Knudsen number 
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E(x, 0) = 2.5 + aw sin —r— 
±j 

with 

a e — 0.3 a u — 0.1 aw = 1- 

This test problem gives rise to a periodic smooth solution in the interval t S [0, 5 x 10 -2 ]. The 
results of this test in log- log scale are shown in Figured] On the left, we reported the stochastic 
error for the pure Monte Carlo and on the right for the Moment Guided method. From top to 
bottom, the errors for the three macroscopic quantities are depicted for different values of the 
Knudsen number. For the Monte Carlo case, the stochastic error does not substantially change 
with respect to the Knudsen number. At variance, for the Moment Guided method, the errors 
decrease as the Knudsen number diminishes. In particular observe that for the test studied the 
stochastic error is between 8 and 12 times smaller of the error of Monte Carlo method for e = 10 -2 
and e = 10~ 3 while is 10~ 4 smaller when e = 10~ 4 . 

5 Conclusions 

We have extended a new class of hybrid methods, first developed in |11) . which aim at reducing 
the variance in Monte Carlo schemes to the general case of the Boltzmann equation. The key idea 
consists in closing the set of macroscopic moments equations through a coupling with the kinetic 
equation solved by means of Monte Carlo methods. In order to correctly close the macroscopic 
system, we need to drive particle positions and velocities in such a way that the moments given 
by the solution of the Boltzmann equation exactly match the moments given by the solution 
of the macroscopic set of moment equations. The new Moment Guided method is constructed 
in such a way that it avoids the problem of stiff regimes or equivalently of very small Knudsen 
numbers. In particular, the method benefits of the asymptotic preserving property which permits 
to circumvent the small scale resolution capturing the limit solution of the problem. A third 
important achievement consists in progressively reduce the contribution of the kinetic solution 
when the Knudsn number tends to zero and to automatically obtain in the fluid limit an high 
order method for the compressible Euler equations without any stochastic contribution. 

The numerical results show reductions of fluctuations in all regimes compared to DSMC method. 
The reduction becomes stronger as we approach equilibrium, when e = the stochastic error be- 
comes automatically zero. Numerical convergence tests show better performances of the proposed 
method compared to pure Monte Carlo schemes for all the cases studied. In particular, for un- 
steady problems, in which time averaging techniques lose their efficiency the Moment Guided 
method seems very promising. It leads to solutions which contain less fluctuations with respect to 
DSMC simulations at a computational cost which is comparable to the cost of a traditional Monte 
Carlo solver with the addition of the cost of a macroscopic solver for the compressible Euler equa- 
tions. This latter is usually computationally much less expensive than any type of solver applied 
to kinetic equations. 

Currently, we are working on extensions of the present method using higher order closures of the 
moment hierarchy in order to solve a larger set of hydrodynamics equations and thus additionally 
reduce the stochastic errors. We are also working to the extension of the method to plasma physics 
problems such as the Vlasov and the collisional Vlasov models. 
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